#include "Sequence.h"


void LucasMatExp(Integer &a, Integer &b, Integer &c, Integer &d, const Integer &exp, const Integer &p, const Integer &q);
void LucasMatExp(Residue &a, Residue &b, Residue &c, Residue &d, const Integer &exp, const Residue &p, const Residue &q);

const Integer Sequence::Factorial(const unsigned n)
{
	if (n <= 92682)
	{
		Integer product(1);

		if (std::IsOdd(n))
			for (unsigned i = 2; i < n; i += 2)
				product *= (i >> 1) * (i + 1);
		else
			for (unsigned i = 4; i <= n; i += 2)
				product *= (i >> 1) * (i - 1);

		return product << (n >> 1);
	}
	else
	{
		Integer p1(1), p2(1);

		for (unsigned i = 4; i <= 92682; i += 2)
			p1 *= (i >> 1) * (i - 1);

		for (unsigned i = 92683; i <= n; ++i)
			p2 *= i;

		return (p1 * p2) << (n >> 1);
	}
}

const Integer Sequence::nCr(const Integer &n, unsigned r)
{
	Integer product(1);

	for (unsigned i = 0; i < r; ++i)
		product = product * (n - i) / (i + 1);

	return product;
}

const Integer Sequence::nPr(const Integer &n, unsigned r)
{
	Integer product(1);

	for (unsigned i = 0; i < r; ++i)
		product = product * (n - i);

	return product;
}

const Integer Sequence::LucasU(const Integer &n, const Integer &p, const Integer &q)
{
	Integer a, b, c, d;
	LucasMatExp(a, b, c, d, n, p, q);
	return c;
}

const Integer Sequence::LucasV(const Integer &n, const Integer &p, const Integer &q)
{
	Integer a, b, c, d;
	LucasMatExp(a, b, c, d, n, p, q);
	return a * 2 + c * p;
}

const Residue Sequence::LucasU(const Integer &n, const Residue &p, const Residue &q)
{
	Residue a, b, c, d;
	LucasMatExp(a, b, c, d, n, p, q);
	return c;
}

const Residue Sequence::LucasV(const Integer &n, const Residue &p, const Residue &q)
{
	Residue a, b, c, d;
	LucasMatExp(a, b, c, d, n, p, q);
	return a * 2 + c * p;
}

void LucasMatExp(Integer &a, Integer &b, Integer &c, Integer &d, const Integer &exp, const Integer &p, const Integer &q)
{
	a = d = 1;
	b = c = 0;

	if (exp == 0)
		return;

	for (int i = exp.Bits(); --i >= 0;)
	{
		const Integer x = b * c;
		const Integer y = a + d;
		a = Square(a) + x;
		d = Square(d) + x;
		b *= y;
		c *= y;

		if (exp.Bit(i))
		{
			Integer temp(std::move(a));
			a = b;
			b = temp * q + b * p;
			std::swap(temp, c);
			c = d;
			d = temp * q + d * p;
		}
	}
}

void LucasMatExp(Residue &a, Residue &b, Residue &c, Residue &d, const Integer &exp, const Residue &p, const Residue &q)
{
	a = d = 1;
	b = c = 0;

	if (exp == 0)
		return;

	for (int i = exp.Bits(); --i >= 0;)
	{
		const Residue x = b * c;
		const Residue y = a + d;
		a = Square(a) + x;
		d = Square(d) + x;
		b *= y;
		c *= y;

		if (exp.Bit(i))
		{
			Residue temp(std::move(a));
			a = b;
			b = temp * q + b * p;
			std::swap(temp, c);
			c = d;
			d = temp * q + d * p;
		}
	}
}
